      Program NLSSTHarmonic
C***********************************************************************
C     This program solves NLS including surface tension.
C     The parameters are set for the Mod expt.
C***********************************************************************
      Parameter(n=1024)
      Complex*16 A,ii,part
      Real*8 h,pi,s,exptdata,delta,epsilon,k0,omega0,gamma
      Integer j,q,carloc
      Dimension A(0:n-1),exptdata(1:3,1:8500)
      Common A
      pi=4.0d0*datan(1.0d0)
      ii=dcmplx(0.0d0,1.0d0)

      open(unit=1,file='ICST.out',status='old')
      read(1,*) exptdata
      close(unit=1)

      carloc=163
      omega0=exptdata(1,carloc)
      k0=1.52951590372984113d0
      epsilon=2.d0*k0*dsqrt(exptdata(2,carloc)**2+exptdata(3,carloc)**2)
      delta=5.122814387346451d0
      s=omega0*2.0d0*pi*epsilon*exptdata(1,2)**(-1)/2.0d0
      gamma=k0**2*71.70d0/981.0d0
      h=k0*epsilon**2/10.0d0
      q=5500

      Do 12 j=0,n-1
         A(j)=0.0d0
 12   Continue

      Do 15 j=0,n-1
         Do 11 jj=-41,41
            part=exptdata(2,carloc+jj)+ii*exptdata(3,carloc+jj)
            A(j)=A(j)+part*exp(jj*ii*2.0d0*pi*j/(n*1.0d0))
 11      Continue
 15   Continue

      Do 13 j=0,n-1
         A(j)=k0/epsilon*A(j)
 13   Continue

      open(unit=2,file="NLSSTHarmonicAlldata.out")
      open(unit=3,file="NLSSTHarmonicAllCQ.out")

      s=s*1.0d0/pi
      Call sixthorder(q,s,h,k0,epsilon,delta,gamma)

      close(2)
      close(3)

      End



      Subroutine linear(s,timestep,k0,delta,epsilon,gamma)
C***********************************************************************
C     This subroutine solves the linear part of 1D NLS.
C***********************************************************************
      Complex*16 A,ii
      Parameter(n=1024)
      Integer i
      Real*8 l,s,timestep,k0,AR,delta,epsilon,gamma,mu1
      Dimension A(0:n-1),AR(0:1,0:n-1)
      Equivalence(A,AR)
      Common A
      ii=dcmplx(0.0d0,1.0d0)

      mu1=(1.0d0+gamma)*(1.0d0-6.0d0*gamma-3.0d0*gamma**2)
      mu1=mu1/(1.0d0+3.0d0*gamma)**3

      Call FFT(AR(0,0),AR(1,0),2,n,1)

      Do 10 i=0,n/2-1
         l=-1.0d0*i/s
         A(i)=A(i)*exp(-mu1*ii*l*l*2.0d0*timestep)
 10   Continue
      
      Do 20 i=n/2+1,n-1
         l=1.0d0*(n-i)/s
         A(i)=A(i)*exp(-mu1*ii*l*l*2.0d0*timestep)
 20   Continue
      
      Call FFT(AR(0,0),AR(1,0),2,n,-1)
      
      Do 30 i=0,n-1
         A(i)=A(i)/(1.0d0*n)
 30   Continue

      Return
      End


      Subroutine nonlinear(timestep,k0,epsilon,s,gamma)
C***********************************************************************
C      This subroutine solves the nonlinear part of 1D NLS.
C**********************************************************************

      Complex*16 A,ii
      Real*8 timestep,k0,epsilon,s,gamma,mu2
      Integer i
      Parameter(n=1024)
      Dimension A(0:n-1)
      Common A
      ii=dcmplx(0.0d0,1.0d0)

      mu2=(8.0d0+gamma+2.0d0*gamma**2)
      mu2=mu2/(2.0d0+2.0d0*gamma-12.0d0*gamma**2)

      Do 10 i=0,n-1
         A(i)=A(i)*exp(mu2*ii*2.0d0*timestep*(cdabs(A(i)))**2)
 10   Continue

      End

      Subroutine sixthorder(q,s,h,k0,epsilon,delta,gamma)
C***********************************************************************
C     Sixth-order operator splitting.
C***********************************************************************
      Complex*16 A
      Real*8 s,h,k0,epsilon,l2n,avewn,delta,qnum,w3,w2,w1,w0,hd2,AR
      Real*8 gamma
      Integer q,j,n
      Parameter(n=1024)
      Parameter(w3=0.7845136104775600d0)
      Parameter(w2=0.2355732133593570d0)
      Parameter(w1=-1.177679984178870d0)
      Parameter(w0=1.3151863206839060d0)
      Dimension A(0:n-1),AR(0:1,0:n-1)
      Equivalence(A,AR)
      Common A

      hd2=h/2.0d0

      Call FFT(AR(0,0),AR(1,0),2,n,1)
      write(2,11,advance="no")(cdabs(A(42-i))/(n*1.0d0),i=1,42)
      write(2,11) (cdabs(A(n-i))/(n*1.0d0),i=1,41)
 11   format(' ',83f9.6)
      Call FFT(AR(0,0),AR(1,0),2,n,-1)
      Do 30 i=0,n-1
         A(i)=A(i)/(1.0d0*n)
 30   Continue
      Call l2norm(l2n)
      Call awn(s,avewn)
      Call qnorm(s,qnum)

      Write(3,5) l2n,avewn,qnum
      Print *,'L2-norm, AWN, and Q',l2n,avewn,qnum
 5    Format(f22.16,f22.16,f22.16)

      Call linear(s,0.50d0*w3*hd2,k0,delta,epsilon,gamma)

      Do 10 j=1,q-1

         Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w3+w2)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w2+w1)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w1+w0)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w0*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w0+w1)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w1+w2)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*(w2+w3)*hd2,k0,delta,epsilon,gamma)
         Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
         Call linear(s,0.50d0*w3*h,k0,delta,epsilon,gamma)
         Call FFT(AR(0,0),AR(1,0),2,n,1)
         write(2,11,advance="no")(cdabs(A(42-i))/(n*1.0d0),i=1,42)
         write(2,11) (cdabs(A(n-i))/(n*1.0d0),i=1,41)
         Call FFT(AR(0,0),AR(1,0),2,n,-1)
         Do 40 i=0,n-1
            A(i)=A(i)/(1.0d0*n)
 40      Continue
         Call l2norm(l2n)
         Call awn(s,avewn)
         Call qnorm(s,qnum)
         Write(3,5) l2n,avewn,qnum
         
 10   Continue

      Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w3+w2)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w2+w1)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w1+w0)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w0*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w0+w1)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w1*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w1+w2)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w2*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*(w2+w3)*hd2,k0,delta,epsilon,gamma)
      Call nonlinear(0.50d0*w3*h,k0,epsilon,s,gamma)
      Call linear(s,0.50d0*w3*hd2,k0,delta,epsilon,gamma)

      Call FFT(AR(0,0),AR(1,0),2,n,1)
      write(2,11,advance="no")(cdabs(A(42-i))/(n*1.0d0),i=1,42)
      write(2,11) (cdabs(A(n-i))/(n*1.0d0),i=1,41)
      Call FFT(AR(0,0),AR(1,0),2,n,-1)
      Do 60 i=0,n-1
         A(i)=A(i)/(1.0d0*n)
 60   Continue

      Call l2norm(l2n)
      Call awn(s,avewn)
      Call qnorm(s,qnum)
      Write(3,5) l2n,avewn,qnum
      Print *,'L2-norm, AWN, and Q',l2n,avewn,qnum

      End



      Subroutine l2norm(l2n)
C***********************************************************************
C     This subroutine calculates the M norm.
C***********************************************************************
      Complex*16 A
      Real*8 l2n
      Integer n,i
      Parameter (n=1024)
      Dimension A(0:n-1)
      Common A

      l2n=0.0d0
      Do 10 i=0,n-1
         l2n=l2n+(cdabs(A(i)))**2
 10   Continue

      l2n=l2n/(n*1.0d0)

      Return
      End

      Subroutine awn(s,avewn)
C***********************************************************************
C     This subroutine calculates the P norm.
C***********************************************************************
      Complex*16 A,B,ii,awnn
      Real*8 pi,s,BR,avewn
      Integer n,x,i
      Parameter(n=1024)
      Dimension A(0:n-1),B(0:n-1),BR(0:1,0:n-1)
      Common A
      Equivalence(B,BR)
      pi=4.0d0*datan(1.0d0)
      ii=dcmplx(0.0d0,1.0d0)

      Do 10 x=0,n-1
         B(x)=A(x)
 10   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,1)
      
      Do 20 i=0,n/2
         B(i)=-ii*i/s*B(i)
 20   Continue
      Do 30 i=n/2+1,n-1
         B(i)=ii*(n-i)/s*B(i)
 30   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,-1)

      Do 40 x=0,n-1
         B(x)=B(x)/(n*1.0d0)
 40   Continue

      awnn=0.0d0

      Do 50 x=0,n-1
         awnn=awnn+A(x)*dconjg(B(x))-B(x)*dconjg(A(x))
 50   Continue

      awnn=0.50d0*ii*awnn/(n*1.0d0)

      avewn=dreal(awnn)

      Return
      End

      Subroutine qnorm(s,qnum)
C***********************************************************************
C     This subroutine calculates the Q norm.
C***********************************************************************
      Complex*16 A,B,ii,qnumm
      Real*8 pi,s,BR,qnum
      Integer n,x,i
      Parameter(n=1024)
      Dimension A(0:n-1),B(0:n-1),BR(0:1,0:n-1)
      Common A
      Equivalence(B,BR)
      pi=4.0d0*datan(1.0d0)
      ii=dcmplx(0.0d0,1.0d0)

      Do 10 x=0,n-1
         B(x)=A(x)
 10   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,1)
      
      Do 20 i=0,n/2
         B(i)=-ii*i/s*B(i)
 20   Continue
      Do 30 i=n/2+1,n-1
         B(i)=ii*(n-i)/s*B(i)
 30   Continue

      Call FFT(BR(0,0),BR(1,0),2,n,-1)

      Do 40 x=0,n-1
         B(x)=B(x)/(n*1.0d0)
 40   Continue

      qnumm=0.0d0
      Do 50 x=0,n-1
         qnumm=qnumm+B(x)*dconjg(B(x))
 50   Continue
      qnum=dreal(qnumm/(n*1.0d0))

      Return
      End

C-----------------------------------------------------------------------
      SUBROUTINE FFT (A,B,IS,N,ID)
C-- +------------------------------------------------------------------+
C-- | A CALL TO  FFT  REPLACES THE COMPLEX DATA VALUES  A(J) + i B(J), |
C-- | J=0,1,..,N-1 WHITH THEIR TRANSFORM                               |
C-- |                                                                  |
C-- |                             2 i ID PI K J / N                    |
C-- |    SUM     (A(K) + i B(K)) e                   ,  J=0,1,..,N-1   |
C-- |  K=0..N-1                                                        |
C-- |                                                                  |
C-- | INPUT AND OUTPUT PARAMETERS                                      |
C-- |     A   ARRAY A(0:*), REAL PART OF DATA/TRANSFORM                |
C-- |     B   ARRAY B(0:*), IMAGINARY PART OF DATA/TRANSFORM           |
C-- | INPUT PARAMETERS                                                 |
C-- |    IS   SPACING IN MEMORY BETWEEN CONSECUTIVE ELEMNTS IN A AND B.|
C---|         (USE IS=+1 FOR ELEMENTS STORED CONSECUTIVELY)            |
C-- |     N   NUMBER OF DATA VALUES, MUST BE A POWER OF TWO            |
C-- |    ID   USE  +1 OR -1  FOR DIRECT AND INVERSE TRANSFORM RESP.    |
C-- |                                                                  |
C-- +------------------------------------------------------------------+
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(0:*),B(0:*)
      J = 0
C---   APPLY PERMUTATION MATRIX   --------
      DO 20 I=0,(N-2)*IS,IS
         IF (I.LT.J) THEN
            TR   = A(J)
            A(J) = A(I)
            A(I) = TR
            TI   = B(J)
            B(J) = B(I)
            B(I) = TI
         ENDIF
         K = IS*N/2
   10    IF (K.LE.J) THEN
            J = J-K
            K = K/2
            GOTO 10
         ENDIF
   20    J = J+K
C---   PERFORM THE  LOG2 N  MATRIX-VECTOR MULT.  ---
      S =  0.0D0
      C = -1.0D0
      L = IS
   30    LH = L
         L  = L+L
         UR = 1.0D0
         UI = 0.0D0
         DO 50 J=0,LH-IS,IS
            DO 40 I=J,(N-1)*IS,L
               IP = I+LH
               TR = A(IP)*UR-B(IP)*UI
               TI = A(IP)*UI+B(IP)*UR
               A(IP) = A(I)-TR
               B(IP) = B(I)-TI
               A(I) = A(I)+TR
   40          B(I) = B(I)+TI
            TI = UR*S+UI*C
            UR = UR*C-UI*S
   50       UI = TI
         S = DSQRT (0.5D0*(1.0D0-C))*ID
         C = DSQRT (0.5D0*(1.0D0+C))
         IF (L.LT.N*IS) GOTO 30
      RETURN
      END
